gusucode.com > 支持向量机工具箱 - LIBSVM OSU_SVM LS_SVM源码程序 > 支持向量机工具箱 - LIBSVM OSU_SVM LS_SVM\stprtool\svm\smomat.m
function [Alpha,bias,nsv,flps]=smomat(aX,aI,aker,aarg,aC,aeps,atol,aAlpha,abias)
% SMOMAT Sequential Minimal Optimization for SVM (L1).
% [Alpha,bias,nsv,flps]=smomat(X,I,ker,arg,C,eps,tol,Alpha, bias)
%
% SMOMAT learns Support Vector Machines (L1) classifier using
% Platt's Sequential Minimal Optimization algorithm.
% It solves binary classification problem with linear
% penalization of classification violations.
%
% The much more faster C-source code version of this
% function is smo.c.
%
% Mandatory inputs:
% X [NxL] L training patterns in N-dimensional space.
% I [1xL] labels of training patterns (1 - 1st, 2 - 2nd class ).
% ker [string] identifier of kernel, see help kernel.
% arg [...] arguments of given kernel, see help kernel.
% C [real] trade-off between margin and training error.
%
% Optional inputs:
% eps [real] tolerance of KKT-conditions fulfilment (default 0.001).
% tol [real] minimal change of optimized Lagrangeians (default 0.001).
% Alpha [1xL] initial values of optimized Lagrangeians. If not given
% then SMO starts from Alpha = zeros(1,L).
% bias [real] initial value of bias. If not given then SMO starts
% from bias = 0.
%
% Outputs:
% Alpha [1xL] found Lagrangeians - solution of the SVM problem.
% bias [real] found threshold of decision function (<w,x>-bias).
% nsv [int] number of support vectors, i.e. number of Alpha > 0.
% flps [int] number of used floating point operatins.
%
% See also SMO, SVMCLASS, SVM.
%
% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
% Written Vojtech Franc (diploma thesis) 23.12.1999, 5.4.2000
% Modifications
% 21-Oct-2001, V.Franc
% 19-September-2001, V. Franc, comments changed.
% 16-September-2001, V. Franc, created
flops(0);
global X Y Alpha b N error_cache tol eps ker arg C;
Y = itosgn(aI); % transform labels (1,2) onto (1,-1)
N = size(aX,2); % number of traning patterns
error_cache = zeros(1,N); % initn error cache
% process the input variables
if nargin < 9,
abias = 0; % threshold
end
if nargin < 8, % use given Lagrangeians or set them to zeros
aAlpha = zeros(1,N);
end
if nargin < 7, % tolerance of KKT-conditions fulfilment for each
atol = 0.001; % pattern
end
if nargin < 6, % minimal change in lagrangeian
aeps = 0.001;
end
% set global variables and clear the locals
X=aX;
b = abias;
Alpha=aAlpha;
tol = atol;
eps = aeps;
arg = aarg;
ker = aker;
C = aC;
clear aX aAlpha atol aeps aarg aker aC;
%-------------------------------------------------------
% MAIN SMO CYCLE
%-------------------------------------------------------
numChanged = 0;
examineAll = 1;
while( numChanged > 0 | examineAll ~= 0 ),
numChanged = 0;
if( examineAll ~= 0 ),
for k=1:N,
numChanged = numChanged + examineExample(k);
end
else
for k=1:N,
if( Alpha(k) ~= 0 & Alpha(k) ~= 0 ),
numChanged = numChanged + examineExample(k);
end
end
end
if( examineAll == 1),
examineAll = 0;
elseif( numChanged == 0 ),
examineAll = 1;
end
end % while( ...)
% we use the following hyperplane <w,x>-bias=0 insted of <w,x>+bias=0
bias = -b;
% clear used global variables
clear global X Y b N error_cache tol eps ker arg C;
nsv = length(find(Alpha>0)); % number of support vectors
flps = flops; % get total number of flops
return; % end of SMO main cycle
% -----------------------------------------------------------
% EXAMINESTEP function
%------------------------------------------------------------
function result=examineExample(i1)
global X Y Alpha b N error_cache tol eps ker arg C;
y1 = Y(i1);
alpha1 = Alpha(i1);
if( alpha1 > 0 & alpha1 < C),
E1 = error_cache(i1); E1 = error_cache(i1);
else
E1 = learned_func(i1) - y1;
end
r1 = y1 * E1;
if((r1 < -tol & alpha1 < C) | (r1 > tol & alpha1 > 0)),
i2 = -1; tmax = 0;
for k=1:N,
if( Alpha(k) > 0 & Alpha(k) < C),
E2 = error_cache(k);
temp = abs(E1 - E2);
if( temp > tmax ),
tmax = temp;
i2 = k;
end
end
end
if( i2 >= 0 ),
if( takeStep (i1, i2) ~= 0),
result = 1;
return;
end
end
randShift = round(N*rand(1));
for k=1:N,
i2 = mod(k-1+randShift,N)+1; % i2 = k;
if( Alpha(i2) > 0 & Alpha(i2) < C),
if(takeStep(i1,i2) ~= 0 ),
result = 1;
return;
end
end
end
randShift = round(N*rand(1));
for k=1:N,
i2 = mod(k-1+randShift,N)+1; % i2 = k;
if(takeStep(i1,i2) ~= 0 ),
result = 1;
return;
end
end
end % if(...
result=0;
return;
% -----------------------------------------------------------
% TAKESTEP function
% -----------------------------------------------------------
function result = takeStep(i1,i2)
global X Y Alpha b N error_cache tol eps ker arg C;
if(i1 == i2),
result = 0;
return;
end
alpha1 = Alpha(i1);
y1 = Y(i1);
if( alpha1 > 0 & alpha1 < C),
E1 = error_cache(i1);
else
E1 = learned_func(i1) - y1;
end
alpha2 = Alpha(i2);
y2 = Y(i2);
if( alpha2 > 0 & alpha2 < C),
E2 = error_cache(i2);
else
E2 = learned_func(i2) - y2;
end
s = y1 * y2;
if(y1 == y2),
gamma = alpha1 + alpha2;
if(gamma > C),
L = gamma-C;
H = C;
else
L = 0;
H = gamma;
end
else
gamma = alpha1 - alpha2;
if(gamma > 0),
L = 0;
H = C - gamma;
else
L = -gamma;
H = C;
end
end
if(L == H),
result = 0;
return;
end
k11 = kernel( X(:,i1), X(:,i1), ker,arg);
k12 = kernel( X(:,i1), X(:,i2), ker,arg);
k22 = kernel( X(:,i2), X(:,i2), ker,arg);
eta = 2 * k12 - k11 - k22;
if( eta < 0 ),
a2 = alpha2 + y2 * (E2 - E1) / eta;
if( a2 < L ),
a2 = L;
elseif( a2 > H ),
a2 = H;
end
else
c1 = eta/2;
c2 = y2 * (E1-E2)- eta * alpha2;
Lobj = c1 * L * L + c2 * L;
Hobj = c1 * H * H + c2 * H;
if( Lobj > Hobj+eps ),
a2 = L;
elseif( Lobj < Hobj-eps ),
a2 = H;
else
a2 = alph2;
end
end
if( a2 < 1e-8),
a2 = 0;
elseif( a2 > C- 1e-8),
a2 = C;
end
if( abs( a2-alpha2 ) < eps*(a2+alpha2+eps)),
result = 0;
return;
end
a1 = alpha1 - s * (a2 - alpha2);
if( a1 < 0 ),
a2 = a2 + s * a1;
a1 = 0;
elseif( a1 > C ),
t = a1-C;
a2 = a2 + s * t;
a1 = C;
end
if( a1 > 0 & a1 < C ),
bnew = b + E1 + y1 * (a1 - alpha1) * k11 + y2 * (a2 - alpha2) * k12;
elseif( a2 > 0 & a2 < C),
bnew = b + E2 + y1 * (a1 - alpha1) * k12 + y2 * (a2 - alpha2) * k22;
else
b1 = b + E1 + y1 * (a1 - alpha1) * k11 + y2 * (a2 - alpha2) * k12;
b2 = b + E2 + y1 * (a1 - alpha1) * k12 + y2 * (a2 - alpha2) * k22;
bnew = (b1 + b2) / 2;
end
delta_b = bnew - b;
b = bnew;
t1 = y1 * (a1-alpha1);
t2 = y2 * (a2-alpha2);
for i=1:N,
if( 0 < Alpha(i) & Alpha(i) < C ),
error_cache(i) = error_cache(i) + t1 * kernel(X(:,i1),X(:,i),ker,arg) + ...
t2 * kernel(X(:,i2),X(:,i),ker,arg) - delta_b;
error_cache(i1) = 0.;
error_cache(i2) = 0.;
end
end
Alpha(i1) = a1;
Alpha(i2) = a2;
result = 1;
return;
%-----------------------------------------------------------
% Compute a value of the learned decision function for
% given patetrn of index k
%-----------------------------------------------------------
function result = learned_func(k)
global X Y Alpha b N error_cache tol eps ker arg C;
result = 0.;
for i=1:N,
if( Alpha(i) > 0),
result = result + Alpha(i)*Y(i)*kernel(X(:,i),X(:,k),ker,arg);
end
end
result = result - b;
return;
% EOF